COVID-19(SARS-CoV-2):

La COVID-19 es la enfermedad causada por el nuevo coronavirus conocido como SARS-CoV-2,el Peru es uno de los paises mas afcetados por esta enfermadad pòr lo que la necesidad de informacion es necesaria par mantener un estricto monitoreo sobre el desarrolo y avance de esta enfermedad y dar a conocer un reflejo cuantitativo de las consecuencias y variables a lo largo de la pandemia.

VARIABLES:

Para el presente trabajo se identificó varias variables a través de la descarga de datos en una serie de tiempo donde presenta una relación directa o indirecta con el desarrollo de la primera y segunda ola del SARS-CoV-2, estas varioables han sido puesta a disposicion de la comunidad y de la cual se descargo los datos de fuentes confiables para desarrrollar el proceso de clusterización.

Esta variables son:

  • Elevación Promedio
  • Índice de Desarrollo Humano (IDH)
  • Varones fallecidos por millón de habitantes
  • Mujeres fallecidas por millón de habitantes
  • Varones contagiados por millón de habitantes
  • Mujeres contagiados por millón de habitantes
  • Adultos por millón de habitantes de habitantes
  • Adultos mayores por millón de habitantes
  • Pobreza por millón de habitantess
  • Pobreza extrema por millón de habitantes
  • Camas UCI por cantidad de hospitales
  • Camas UCIN por cantidad de hospitales
  • Ventiladores mecánicos por cantidad de hospitales
library(dplyr)
library(tidyverse)
library(ggplot2)
library(car)
library(dplyr)
library(factoextra)
library(cluster)
library(simstudy)
library(data.table)
library(PerformanceAnalytics)
library(ggdark)
library(scatterplot3d)
library( rgl )
library(abind)
library(magic)
library(igraph)
library(LICORS)
library(clValid)
library(clusterSim)
library(ggthemes)
library(ggthemes)
library(sf)
dataset1 <- read_csv("https://raw.githubusercontent.com/ElizabethRamosCastillo/PROGRAMACION_R/master/PRESENTACIONES/PROYECTO_FINAL/data/Data_final/Dataset_1_Ola.csv")
## Rows: 196 Columns: 14
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (1): PROVINCIA
## dbl (13): ElevacionPromedio, IDH_2019, FallMasc_1ola_mill_hab, FallFem_1ola_...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
dataset2 <- read_csv("https://raw.githubusercontent.com/ElizabethRamosCastillo/PROGRAMACION_R/master/PRESENTACIONES/PROYECTO_FINAL/data/Data_final/Dataset_2_Ola.csv")
## Rows: 196 Columns: 14
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr  (1): PROVINCIA
## dbl (13): ElevacionPromedio, IDH_2019, FallMasc_2ola_mill_hab, FallFem_2ola_...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

##PRIMERA OLA

HISTOGRAMAS

  1. Histograma de Elevación Promedio:
hist_EP <- 
  ggplot(data = dataset1,
         mapping = aes(x = ElevacionPromedio)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de Elevación Promedio") + 
  labs(x ="Elevación Promedio") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
## Inverted geom defaults of fill and color/colour.
## To change them back, use invert_geom_defaults().
hist_EP

  1. Histograma de Índice de Desarrollo Humano (IDH):
hist_IDH <- 
  ggplot(data = dataset1,
         mapping = aes(x = IDH_2019)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de Índice de Desarrollo Humano") + 
  labs(x ="IDH") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_IDH 

#3. Histograma de Varones fallecidos por millón de habitantes:

hist_Var_fall <- 
  ggplot(data = dataset1,
         mapping = aes(x = FallMasc_1ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de varones fallecidos por millón de habitantes") + 
  labs(x ="Varones fallecidos por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Var_fall

#4. Histograma de Mujeres fallecidas por millón de habitantes:

hist_Muj_fall <- 
  ggplot(data = dataset1,
         mapping = aes(x = FallFem_1ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de mujeres fallecidas por por millón de habitantes") + 
  labs(x ="Mujeres fallecidos por por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Muj_fall

#5. Histograma de Varones contagiados por millón de habitantes:

hist_Var_cont <- 
  ggplot(data = dataset1,
         mapping = aes(x = CasosMasc_1ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de varones contagiados por millón de habitantes") + 
  labs(x ="Varones contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Var_cont 

#6. Histograma de mujeres contagiados por millón de habitantes:

hist_Muj_cont <- 
  ggplot(data = dataset1,
         mapping = aes(x = CasosFem_1ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de mujeres contagiados por millón de habitantes") + 
  labs(x ="Mujeres contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Muj_cont

#7. Histograma de adultos por millón de habitantes:

hist_Adult <- 
  ggplot(data = dataset1,
         mapping = aes(x = Adult_15_64_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de adultos por millón de habitantes") + 
  labs(x ="Adultos contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Adult

#8. Histograma de adultos mayores por millón de habitantes:

hist_Adult_mayor <- 
  ggplot(data = dataset1,
         mapping = aes(x = Mayores_65_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de adultos mayores por millón de habitantes") + 
  labs(x ="Adultos mayores contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Adult_mayor 

#9. Histograma de pobreza por millón de habitantes:

hist_Pobrez <- 
  ggplot(data = dataset1,
         mapping = aes(x = Pobreza_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de pobreza por millón de habitantes") + 
  labs(x ="Pobreza por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Pobrez

#10. Histograma de pobreza extrema por millón de habitantes:

hist_Pob_ext <- 
  ggplot(data = dataset1,
  mapping = aes(x = PobrezaExt_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de pobreza extrema por millón de habitantes") + 
  labs(x ="Pobreza extrema por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Pob_ext

#11. Histograma de camas UCI por cantidad de hospitales:

hist_UCI <- 
  ggplot(data = dataset1,
         mapping = aes(x = UCI_1ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de camas UCI por cantidad de hospitales") + 
  labs(x ="Camas UCI por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_UCI

#12. Histograma de camas UCIN por cantidad de hospitales:

hist_UCIN <- 
  ggplot(data = dataset1,
         mapping = aes(x = UCIN_1ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de camas UCIN por cantidad de hospitales") + 
  labs(x ="Camas UCIN por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_UCIN 

#13. Histograma de ventiladores mecánicos por cantidad de hospitales:

hist_VM <- 
  ggplot(data = dataset1,
         mapping = aes(x = VENT_1ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de ventiladores mecánicos por cantidad de hospitales") + 
  labs(x ="Ventiladores Mecánicos por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_VM

SEGUNDA OLA

HISTOGRAMAS

library(tidyverse)
library(ggdark)
library(scatterplot3d)
library( rgl )
library(abind)
library(magic)
  1. Histograma de Elevación Promedio:
hist_EP <- 
  ggplot(data = dataset2,
         mapping = aes(x = ElevacionPromedio)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de Elevación promedio") + 
  labs(x ="Elevación Promedio") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_EP

#2. Histograma de Índice de Desarrollo Humano (IDH):

hist_IDH <- 
  ggplot(data = dataset2,
         mapping = aes(x = IDH_2019)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de Índice de Desarrollo Humano") + 
  labs(x ="IDH") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_IDH

#3. Histograma de Varones fallecidos por millón de habitantes:

hist_Var_fall <- 
  ggplot(data = dataset2,
         mapping = aes(x = FallMasc_2ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de varones fallecidos por millón de habitantes") + 
  labs(x ="Varones fallecidos por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Var_fall

#4. Histograma de Mujeres fallecidas por millón de habitantes:

hist_Muj_fall <- 
  ggplot(data = dataset2,
         mapping = aes(x = FallFem_2ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de mujeres fallecidas por millón de habitantes") + 
  labs(x ="Mujeres fallecidos por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Muj_fall

#5. Histograma de Varones contagiados por millón de habitantes:

hist_Var_cont <- 
  ggplot(data = dataset2,
         mapping = aes(x = CasosMasc_2ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de varones contagiados por millón de habitantes") + 
  labs(x ="Varones contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Var_cont

#6. Histograma de mujeres contagiados por millón de habitantes:

hist_Muj_cont <- 
  ggplot(data = dataset2,
         mapping = aes(x = CasosFem_2ola_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de mujeres contagiados por millón de habitantes") + 
  labs(x ="Mujeres contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Muj_cont

#7. Histograma de adultos por millón de habitantes:

hist_Adult <- 
  ggplot(data = dataset2,
         mapping = aes(x = Adult_15_64_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de adultos por millón de habitantes") + 
  labs(x ="Adultos contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Adult

#8. Histograma de adultos mayores por millón de habitantes:

hist_Adult_mayor <- 
  ggplot(data = dataset2,
         mapping = aes(x = Mayores_65_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de adultos mayores por millón de habitantes") + 
  labs(x ="Adultos mayores contagiados por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Adult_mayor

#9. Histograma de pobreza por millón de habitantes:

hist_Pobrez <- 
  ggplot(data = dataset2,
         mapping = aes(x = Pobreza_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de pobreza por millón de habitantes") + 
  labs(x ="Pobreza por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Pobrez

#10. Histograma de pobreza extrema por millón de habitantess:

hist_Pob_ext <- 
  ggplot(data = dataset2,
         mapping = aes(x = PobrezaExt_mill_hab)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de pobreza extrema por millón de habitantes") + 
  labs(x ="Pobreza extrema por millón de habitantes") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_Pob_ext

#11. Histograma de camas UCI por cantidad de hospitales:

hist_UCI <- 
  ggplot(data = dataset2,
         mapping = aes(x = UCI_2ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de camas UCI por cantidad de hospitales") + 
  labs(x ="Camas UCI por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_UCI 

#12. Histograma de camas UCIN por cantidad de hospitales:

hist_UCIN <- 
  ggplot(data = dataset2,
         mapping = aes(x = UCIN_2ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de camas UCIN por cantidad de hospitales") + 
  labs(x ="Camas UCIN por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_UCIN

#13. Histograma de ventiladores mecánicos por cantidad de hospitales:

hist_VM <- 
  ggplot(data = dataset2,
         mapping = aes(x = VENT_2ola_cant_hosp)) +
  geom_histogram(bins = 9, color="black", fill = "orange" ) +
  dark_theme_gray()+
  ggtitle("Histograma de ventiladores mecánicos por cantidad de hospitales") + 
  labs(x ="Ventiladores Mecánicos por hospitales") +
  theme(
    plot.title = 
      element_text(vjust = 1, 
                   hjust = 0.5,
                   size=rel(1.7),
                   face="bold",
                   color = "red")) +
  theme(axis.title.x = element_text(face="bold", vjust=-0.5, colour="white", size=rel(1.5))) +
  theme(axis.title.y = element_text(face="bold", vjust=1.5, colour="white", size=rel(1.5)))
hist_VM

DESARROLLO:

Removemos:

rm(list = ls())

Cargamos las librerias con pacman, esto no te notifica los warnings, library si te sale.

pacman::p_load(tidyverse) View(Dataset_1_Ola)

Cargamos la data:

  • Lectura de archivos
Data_1_Ola <- read.csv("https://raw.githubusercontent.com/ElizabethRamosCastillo/PROGRAMACION_R/master/PRESENTACIONES/PROYECTO_FINAL/data/Data_final/Dataset_1_Ola.csv")
rownames(Data_1_Ola) = Data_1_Ola$PROVINCIA

View(Data_1_Ola)

ANÁLISIS EXPLORATORIO DE DATOS (EDA)

VAMOS A IR GENERANDO PRIMERO GRAFICOS UNIVARIADOS:

ELEVACION

library(dplyr)
library(tidyverse)
library(ggplot2)
library(car)
library(dplyr)
library(factoextra)
library(cluster)
library(simstudy)
library(data.table)
library(PerformanceAnalytics)
plot(
  Data_1_Ola$ElevacionPromedio, 
  pch = 20
)

text( #Para colocar el texto.
  Data_1_Ola$ElevacionPromedio,
  labels = Data_1_Ola$PROVINCIA,
  cex = 0.5, col = "blue", pos = 1
)

CASOS DE COVID EN VARONES

plot(
  Data_1_Ola$FallMasc_1ola_mill_hab, 
  pch = 20
)

text( #Para colocar el texto.
  Data_1_Ola$FallMasc_1ola_mill_hab,
  labels = Data_1_Ola$PROVINCIA,
  cex = 0.5, col = "blue", pos = 1
)

GRÁFICOS BIVARIADOS:

ELEVACIÓN - FALLECIDOS MASCULINOS

plot(
  Data_1_Ola$ElevacionPromedio, 
  Data_1_Ola$FallMasc_1ola_mill_hab,
  pch = 20, col = "red"
)

text( #Para colocar el texto.
  Data_1_Ola$ElevacionPromedio, 
  Data_1_Ola$FallMasc_1ola_mill_hab,
  labels = Data_1_Ola$PROVINCIA,
  cex = 0.5, col = "blue", pos = 1
)

GRÁFICO TRIDIMENSIONAL:

Se utilizará la librería scatterplot3d

library(scatterplot3d)
s3d<- scatterplot3d(
  Data_1_Ola$ElevacionPromedio, 
  Data_1_Ola$FallMasc_1ola_mill_hab,
  Data_1_Ola$FallFem_1ola_mill_hab,
  angle = 40, 
  color = "blue",
  pch = 20
)

text(
  s3d$xyz.convert(
    Data_1_Ola$ElevacionPromedio, 
    Data_1_Ola$FallMasc_1ola_mill_hab,
    Data_1_Ola$FallFem_1ola_mill_hab
  ),
  labels = Data_1_Ola$PROVINCIA,
  cex = 0.5,
  pos = 4
)

ESTANDARIZACIÓN DE LOS DATOS:

NOSOTROS TENEMOS QUE ESTARANDIZAR NUESTROS DATOS, ES UNA PARTE PREVIO A LOS PCA, PORQUE TENEMOS VARIAS VARIABLES.

UNA VARIABLE ESTANDARIZADA TIENE LA PROPIEDAD DE QUE LA DIAGONAL ES 1, SE HACE CON LA FUNCION SCALE PARA ESTANDARIZAR

Data_1_Ola_PCA <- Data_1_Ola %>% 
  dplyr::select(-PROVINCIA)

Estandarizamos los datos: (Normalización)

scale_data_1ola <- scale(Data_1_Ola_PCA) #Te estandariza los datos
View(scale_data_1ola)

Calculando la matriz varianza covarianza

cov_1ola<-cov(scale_data_1ola) #La diagonal representa a la varianza de cada variable y los demas valores representan a las covarianza.
View(cov_1ola)

La suma de los elementos de la diagonal es 7 (que es la varianza explicada global de los datos), osea la suma de autovalores sera igual a 7 La media de todas las varianzas seran 1.

diag(cov_1ola) %>% sum()
## [1] 13

La suma de las varianzas son conocido como la TRAZA DE LA MATRIZ, ENTONCES LA SUMA DE LA VARIANZAS SERA IGUAL AL NUMERO DE VARIABLES.

la función describe te saca los estadásticos de las variables:

library(psych)

describe(Data_1_Ola_PCA)

JUSTIFICACION CON UNA MATRIZ DE CORRELACIONES

library(PerformanceAnalytics)#para sacar la correlacion
library(GGally)
library(lattice)
library(cluster)
library(ggplot2)
library(corrr)
library(cptcity)#para colores

Primero calculamos la matriz de correlacion de los “datos estandarizados”

mtx_1ola <- cor(scale_data_1ola)#introducimos el dataset estandarizado
chart.Correlation(scale_data_1ola, histogram = T, pch = 20)

Los asteristicos significan el nivel de significacion.

ggpairs(as.data.frame(scale_data_1ola))#solo te aceptan data.frame, no matrices 

library(corrplot)
corrplot(mtx_1ola, method ="shade", sade.col=NA, tl.col ="black", tl.srt =35,
          addCoef.col = "black", addcolorlabel ="no", order ="AOE",tl.cex = 0.5)
## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "sade.col" is not a graphical parameter
## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "addcolorlabel" is not a graphical parameter
## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "sade.col" is not a graphical parameter
## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "addcolorlabel" is not a graphical parameter
## Warning in title(title, ...): "sade.col" is not a graphical parameter
## Warning in title(title, ...): "addcolorlabel" is not a graphical parameter

ANALISIS DE COMPONENTES PRINCIPALES:

Cargar el paquete ade4 para realizar el analisis de componentes principales

library(ade4)

El scale te dice si quieres estandaraizar los datos, en este caso sera false porque ya lo estanfarizamos El scanmf te explica que si deseas plotear en grafico de SEDIMENTACION. nf es el numero de componentes o ejes retenido, esto es igual a la cantidad de variables como maximo, pero se puede retener menos

PCA_1ola <-dudi.pca(scale_data_1ola , #la funcion que hace PCA es dudi.pca
               scale = F, scannf = F, 
               nf = ncol(scale_data_1ola))
summary(PCA_1ola)
## Class: pca dudi
## Call: dudi.pca(df = scale_data_1ola, scale = F, scannf = F, nf = ncol(scale_data_1ola))
## 
## Total inertia: 12.93
## 
## Eigenvalues:
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  6.3727  1.7135  1.5145  1.0406  0.7721 
## 
## Projected inertia (%):
##     Ax1     Ax2     Ax3     Ax4     Ax5 
##  49.272  13.248  11.710   8.046   5.970 
## 
## Cumulative projected inertia (%):
##     Ax1   Ax1:2   Ax1:3   Ax1:4   Ax1:5 
##   49.27   62.52   74.23   82.28   88.25 
## 
## (Only 5 dimensions (out of 13) are shown)

EIGENVALUES: Nos devuelve los autovalores de los componentes (es como un head que te muestra 5 primeros)

PROJECT INTERTIA: Luego te presenta la varianza explicada CUMULATIVE PROJECTED INTERTIA: Luego te presenta la varianza explicada acumulada.

PRIMER CRITERIO: AUTOVALORES

PARA OBTENER LOS VALORES PROPIOS (AUTOVALORES)

PCA_1ola$eig
##  [1] 6.37272527 1.71351508 1.51452516 1.04058940 0.77209142 0.50150165
##  [7] 0.31517357 0.24632625 0.19512874 0.11080086 0.09266958 0.04275247
## [13] 0.01587401
sum(PCA_1ola$eig) #CUANDO SE SUMA LOS AUTOVALORES TE DA EL NUMERO DE VARIABLES QUE ES 13, PERO COMO NADA ES PRECISO, TE DA UNA APROXIMACIÓN
## [1] 12.93367

CON ESTOS AUTOVALORES PUEDES DETERMINAR EL PRIMER CRITERIO PARA OBTENER LOS VALORES VECTORES O AUTOVECTORES

PCA_1ola$c1 #los vectores propios representan los laudings, es decir los "W",  o tambien llamado los pesos,estos varian por variable y por componentes
##                                CS1           CS2         CS3         CS4
## ElevacionPromedio        0.1798862  0.4926498371  0.07858427 -0.02576456
## IDH_2019                -0.3408322  0.1068745493  0.24982652  0.12825503
## FallMasc_1ola_mill_hab  -0.3513863 -0.0846709952  0.08808144 -0.22970915
## FallFem_1ola_mill_hab   -0.2989639  0.0006817441  0.20481255 -0.29355496
## CasosMasc_1ola_mill_hab -0.2973666 -0.3521644490 -0.04943442 -0.17004217
## CasosFem_1ola_mill_hab  -0.2774840 -0.4046544554 -0.08334619 -0.17619683
## Adult_15_64_mill_hab    -0.2776572  0.1092564004  0.09166543 -0.27170083
## Mayores_65_mill_hab      0.0744949  0.4181276765  0.29450059 -0.62398548
## Pobreza_mill_hab         0.2938549 -0.1372391680 -0.37159663 -0.33649459
## PobrezaExt_mill_hab      0.2627519 -0.1131171784 -0.38311625 -0.43771326
## UCI_1ola_cant_hosp      -0.2807647  0.3028717538 -0.42503742  0.07022414
## UCIN_1ola_cant_hosp     -0.2863531  0.2057363339 -0.35965245 -0.04030107
## VENT_1ola_cant_hosp     -0.2726521  0.3138196306 -0.42698056  0.09710836
##                                  CS5         CS6         CS7         CS8
## ElevacionPromedio        0.557841476 -0.02757129  0.59296575  0.20093203
## IDH_2019                -0.007092714 -0.14004957 -0.09238350  0.01376427
## FallMasc_1ola_mill_hab  -0.177945066  0.16782453  0.16804029  0.11051254
## FallFem_1ola_mill_hab   -0.348174856  0.29139629  0.55323448 -0.27248160
## CasosMasc_1ola_mill_hab  0.474333571  0.03908436 -0.04633604  0.06170387
## CasosFem_1ola_mill_hab   0.464327233  0.08123052 -0.01383812  0.01939185
## Adult_15_64_mill_hab    -0.078005632 -0.86890748  0.02524939  0.01295563
## Mayores_65_mill_hab      0.118192428  0.24101340 -0.51049440 -0.03273144
## Pobreza_mill_hab        -0.115881087 -0.09585542  0.10863508  0.08885015
## PobrezaExt_mill_hab     -0.060604652 -0.09726541  0.11178912 -0.14988715
## UCI_1ola_cant_hosp       0.059506758  0.08064326 -0.12348035 -0.30013984
## UCIN_1ola_cant_hosp     -0.215098181  0.14086347 -0.04646875  0.77229156
## VENT_1ola_cant_hosp      0.091193900  0.01845925 -0.01628697 -0.38817049
##                                 CS9        CS10          CS11         CS12
## ElevacionPromedio        0.03001878  0.09601215  0.0007122414  0.058445150
## IDH_2019                 0.62615503  0.19446651 -0.5653798913  0.088032376
## FallMasc_1ola_mill_hab  -0.13108931  0.79054413  0.2540895636  0.016849969
## FallFem_1ola_mill_hab   -0.01667383 -0.41969623 -0.1290532271  0.016151185
## CasosMasc_1ola_mill_hab  0.02837189 -0.05663003 -0.0531355720  0.012754804
## CasosFem_1ola_mill_hab  -0.05662384 -0.15682140 -0.0227890703 -0.028286936
## Adult_15_64_mill_hab    -0.18750545 -0.11612194  0.1343353826  0.004815958
## Mayores_65_mill_hab     -0.05728607 -0.04403264 -0.0400259404 -0.052113774
## Pobreza_mill_hab        -0.29720556  0.18601787 -0.6820785950  0.114048464
## PobrezaExt_mill_hab      0.66308024  0.03829466  0.2990389183 -0.031532744
## UCI_1ola_cant_hosp      -0.06778197 -0.01324996  0.0808916937  0.719050103
## UCIN_1ola_cant_hosp      0.09015113 -0.24905335  0.0385087308 -0.091254194
## VENT_1ola_cant_hosp     -0.07991830  0.10180448 -0.1160006677 -0.667252624
##                                 CS13
## ElevacionPromedio        0.025152455
## IDH_2019                 0.088706122
## FallMasc_1ola_mill_hab   0.039546403
## FallFem_1ola_mill_hab   -0.052026542
## CasosMasc_1ola_mill_hab -0.718834392
## CasosFem_1ola_mill_hab   0.685113908
## Adult_15_64_mill_hab     0.011972935
## Mayores_65_mill_hab      0.016124521
## Pobreza_mill_hab        -0.004169489
## PobrezaExt_mill_hab      0.011750432
## UCI_1ola_cant_hosp       0.011379375
## UCIN_1ola_cant_hosp      0.005656789
## VENT_1ola_cant_hosp     -0.020091057

EXPLICACIÓN DE RETENCION LAS VARIABLES:

EL PRIMER CRITERIO PARA SABER CUANTAS COMPONENTES VAMOS A RETENER SE BASA EN LOS AUTOVALORES: DICE QUE VAMOS A RETENER TANTOS COMPONENTES CONFORME ESTOS TENGAN ASOCIADOS AUTOVALORES MAYORES A 1 PORQUE ESTE CRITERIO ES CONOCIDO COMO LA MEDIA ARIMETICA, CONFORMEN LOS AUTOVALORES SEAN MAYORES AL PROMEDIO DE LA VARIANZA PROMEDIO QUE ERA 1, YA QUE NOSOTROS ESTANDARIZAMOS LAS VARIABLES, Y ESTO HACIA QUE LAS DIAGONALES ERAN 1, Y SU PROMEDIO SERA 1, PERO SI NO FUERON ESTANDARIZADAS, CALCULAS LA MEDIA ARIMETICA Y VALORES MAYORES A LA MEDIA ARIMETICA SERAN LOS QUE RETENGAMOS. ENTONCES NOS QUEDAMOS CON LAS 4 PRIMERAS COMPONENTES

TERCER CRITERIO: RETENCIÓN DE VARIABLES

TERCER CRITERIO: CRITERIO DE RETENCION DE VARIABLES:

SE PUEDEN ELIMINAR VARIABLES, PERO COMO SE ELIMINA, EN FUNCION A LA MATRIZ DE CORRELACION ENTRE LAS VARIABLES ORIGINALES Y LAS COMPONENTES DE LA MATRIZ.

PLOTEANDO LA MATRIZ COMPONENTES VS VARIABLES ORIGINALES

PCA_1ola$co 
##                              Comp1         Comp2       Comp3       Comp4
## ElevacionPromedio        0.4541093  0.6448850505  0.09671055 -0.02628224
## IDH_2019                -0.8604056  0.1399001763  0.30745162  0.13083204
## FallMasc_1ola_mill_hab  -0.8870486 -0.1108354351  0.10839834 -0.23432466
## FallFem_1ola_mill_hab   -0.7547121  0.0008924119  0.25205471 -0.29945332
## CasosMasc_1ola_mill_hab -0.7506799 -0.4609878484 -0.06083699 -0.17345880
## CasosFem_1ola_mill_hab  -0.7004878 -0.5296979502 -0.10257086 -0.17973712
## Adult_15_64_mill_hab    -0.7009251  0.1430180506  0.11280901 -0.27716007
## Mayores_65_mill_hab      0.1880568  0.5473345721  0.36243023 -0.63652312
## Pobreza_mill_hab         0.7418148 -0.1796478576 -0.45730927 -0.34325572
## PobrezaExt_mill_hab      0.6632976 -0.1480718591 -0.47148602 -0.44650816
## UCI_1ola_cant_hosp      -0.7087695  0.3964630687 -0.52307675  0.07163514
## UCIN_1ola_cant_hosp     -0.7228772  0.2693115395 -0.44261004 -0.04111084
## VENT_1ola_cant_hosp     -0.6882901  0.4107939820 -0.52546809  0.09905954
##                                Comp5       Comp6        Comp7        Comp8
## ElevacionPromedio        0.490168235 -0.01952510  0.332892876  0.099725109
## IDH_2019                -0.006232278 -0.09917860 -0.051864394  0.006831382
## FallMasc_1ola_mill_hab  -0.156358074  0.11884793  0.094338361  0.054848770
## FallFem_1ola_mill_hab   -0.305936832  0.20635747  0.310587616 -0.135236067
## CasosMasc_1ola_mill_hab  0.416790897  0.02767829 -0.026013202  0.030624411
## CasosFem_1ola_mill_hab   0.407998454  0.05752484 -0.007768768  0.009624419
## Adult_15_64_mill_hab    -0.068542560 -0.61533231  0.014175091  0.006430045
## Mayores_65_mill_hab      0.103854188  0.17067794 -0.286593195 -0.016245027
## Pobreza_mill_hab        -0.101823243 -0.06788172  0.060988083  0.044097456
## PobrezaExt_mill_hab     -0.053252540 -0.06888023  0.062758771 -0.074390890
## UCI_1ola_cant_hosp       0.052287834  0.05710896 -0.069322263 -0.148963200
## UCIN_1ola_cant_hosp     -0.189004046  0.09975497 -0.026087703  0.383298076
## VENT_1ola_cant_hosp      0.080130924  0.01307225 -0.009143556 -0.192653925
##                                Comp9       Comp10        Comp11        Comp12
## ElevacionPromedio        0.013260308  0.031959336  0.0002168182  0.0120845117
## IDH_2019                 0.276593837  0.064731608 -0.1721111276  0.0182021652
## FallMasc_1ola_mill_hab  -0.057906580  0.263146555  0.0773491275  0.0034840128
## FallFem_1ola_mill_hab   -0.007365393 -0.139703292 -0.0392859682  0.0033395274
## CasosMasc_1ola_mill_hab  0.012532822 -0.018850305 -0.0161753599  0.0026372688
## CasosFem_1ola_mill_hab  -0.025012663 -0.052200768 -0.0069373755 -0.0058487969
## Adult_15_64_mill_hab    -0.082827496 -0.038653234  0.0408939450  0.0009957798
## Mayores_65_mill_hab     -0.025305194 -0.014657039 -0.0121845680 -0.0107753939
## Pobreza_mill_hab        -0.131285742  0.061919328 -0.2076361715  0.0235814264
## PobrezaExt_mill_hab      0.292904952  0.012747052  0.0910324654 -0.0065199219
## UCI_1ola_cant_hosp      -0.029941585 -0.004410482  0.0246247892  0.1486756281
## UCIN_1ola_cant_hosp      0.039822802 -0.082901799  0.0117227039 -0.0188683300
## VENT_1ola_cant_hosp     -0.035302613  0.033887415 -0.0353125501 -0.1379656334
##                                Comp13
## ElevacionPromedio        0.0031690108
## IDH_2019                 0.0111762713
## FallMasc_1ola_mill_hab   0.0049825346
## FallFem_1ola_mill_hab   -0.0065549337
## CasosMasc_1ola_mill_hab -0.0905674600
## CasosFem_1ola_mill_hab   0.0863189452
## Adult_15_64_mill_hab     0.0015084953
## Mayores_65_mill_hab      0.0020315624
## Pobreza_mill_hab        -0.0005253227
## PobrezaExt_mill_hab      0.0014804617
## UCI_1ola_cant_hosp       0.0014337115
## UCIN_1ola_cant_hosp      0.0007127108
## VENT_1ola_cant_hosp     -0.0025313146
PCA_PLOT_1ola<-PCA_1ola$co

corrplot(as.matrix(PCA_PLOT_1ola),method ="shade",addCoef.col = "black", addcolorlabel ="no",tl.cex = 0.9)
## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "addcolorlabel" is not a graphical parameter
## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "addcolorlabel" is not a graphical parameter
## Warning in title(title, ...): "addcolorlabel" is not a graphical parameter

Donde:

cw: los pesos de la columna lw: los pesos de la fila eig: los valores propios rank: el rango de la matriz analizada nf: el número de factores guardados c1: las puntuaciones normalizadas de la columna, es decir, los ejes principales l1:las puntuaciones normativas de la fila co:las coordenadas de la columna li: las coordenadas de la fila, es decir, los componentes principales call: la función de llamada cent: el vector p que contiene las medias de las variables (tenga en cuenta que si centro = F, el vector contiene p 0) norm: el vector p que contiene las desviaciones estándar de las variables, es decir, la raíz de la suma de los cuadrados, las desviaciones de los valores de sus medias divididas por n (tenga en cuenta que si norma = F, el vector contiene p 1)

ENTONCES DE ESTE ANALISIS VISUALIZAMOS SI ALGUNAS DE LAS VARIABLES NO SE RELACIONAN CON ALGUNA COMPONENTE, EN ESTE CASO VEMOS QUE SI SE RELACIONAN TODOS, DONDE X1,X2,…X7 SON LAS VARIABLES ORIGINALES Y COMP1, COMP2, ETC SON LOS COMPONENTES PRINCIPALES.

SI VEMOS QUE UNA VARIABLE NO SE CORRELACIONA POR COMPLETO, LE QUITAMOS ESA VARIABLE DE LA MATRIZ ORIGINAL Y VOLVEMOS A CORRER DENUEVO TODO EL ANALISIS PCA. CONCLUSION: CON ESTA MATRIZ VEMOS QUE LAS CORRELACIONES SON MAYORES DE +/-0.5 EN LA PRIMERA Y SEGUNDA COMPONENTE, POR LO QUE SE RETENDRA LAS 2 PRIMERAS COMPONENTES.

CONTRIBUCION DE VARIABLES A LAS COMPONENTES PRINCIPALES

contrib_1ola <- as.matrix(PCA_1ola$co * PCA_1ola$co)

Si elevamos al cuadrado las correlaciones, tendremos una mejor representacion, un valor alto representa una contribucion alta, Y un resultado bajo representa una contribucion baja

corrplot(contrib_1ola, is.corr = F,method ="shade",addCoef.col = "black", addcolorlabel ="no",tl.cex = 0.9)
## Warning in text.default(pos.xlabel[, 1], pos.xlabel[, 2], newcolnames, srt =
## tl.srt, : "addcolorlabel" is not a graphical parameter
## Warning in text.default(pos.ylabel[, 1], pos.ylabel[, 2], newrownames, col =
## tl.col, : "addcolorlabel" is not a graphical parameter
## Warning in title(title, ...): "addcolorlabel" is not a graphical parameter

is.corr es para que no correlacione desde -1 a +1, sino de 0 a 1, dado que los resultados estan elevados al cuadrado (valores positivos)

#PODEMOS CONCLUIR QUE LAS VARIABLES DE PERSONAS MAYORES A 65 AÑOS Y ELEVACIÓN CONTRIBUYEN MEJOR A LA SEGUNDA COMPONENTE LAS DEMÁS VARIABLES CONTRIBUYEN MEJOR EN LA PRIEMRA COMPONENTE #SEGUN ESTE GRAFICO, SOLO RETENDRIAMOS LAS 2 PRIMERAS COMPONENTES, DADO QUE LA COMPONENTE 1 Y 2 EXPLICAN MEJOR CADA UNA DE LAS 13 VARIABLES QUE LAS DEMÁS COMPONENTES.

OBTENCIÓN DE LOS SCORES O PUNTUACIONES:

LA SALIDA DEL PCA TIENE QUE HACER LO MISMO, UNA TABLA DE 196x13

as.tibble(scale_data_1ola)
## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
## # A tibble: 196 x 13
##    ElevacionPromedio IDH_2019 FallMasc_1ola_m~ FallFem_1ola_mi~ CasosMasc_1ola_~
##                <dbl>    <dbl>            <dbl>            <dbl>            <dbl>
##  1             0.446   0.653            -0.174           0.0579           0.430 
##  2             0.624  -1.16             -0.634          -0.760           -0.640 
##  3             1.07   -1.25             -0.656          -0.295           -0.761 
##  4             0.177  -0.641            -0.779          -1.01            -0.288 
##  5            -1.25   -0.0370            0.389           0.518           -0.0255
##  6             0.624  -0.555            -0.693          -0.494           -0.638 
##  7             0.356  -0.210            -0.787          -0.572           -0.569 
##  8             0.893  -0.986            -0.356          -0.388           -0.394 
##  9             0.995  -0.468            -0.430          -0.130           -0.571 
## 10             1.25   -1.07             -0.818          -0.862           -0.696 
## # ... with 186 more rows, and 8 more variables: CasosFem_1ola_mill_hab <dbl>,
## #   Adult_15_64_mill_hab <dbl>, Mayores_65_mill_hab <dbl>,
## #   Pobreza_mill_hab <dbl>, PobrezaExt_mill_hab <dbl>,
## #   UCI_1ola_cant_hosp <dbl>, UCIN_1ola_cant_hosp <dbl>,
## #   VENT_1ola_cant_hosp <dbl>

ENTONCES VAMOS A TENER LA NUEVA TABLA DE PUNTUACIONES O SCORES

PCA_PRIMERA_OLA <- PCA_1ola$li #estoy valores ya son los valores resultantes del PCA, el cual se les conocen como SCORES o Puntuaciones
head(PCA_PRIMERA_OLA)
##                    Axis1      Axis2      Axis3      Axis4      Axis5      Axis6
## ABANCAY       -3.4919948  2.7915258 -3.1313116  0.2103763  0.7547921 -0.4381602
## ACOBAMBA       2.0963553 -0.4770106 -0.1553921  2.2149076  0.3234739  1.6660686
## ACOMAYO        2.1663012  0.7270715  0.3506939 -0.1628312  0.2269999  0.2316588
## AIJA           1.8371843  0.4596678  0.6078255 -0.1877961  0.4789561  0.3981455
## ALTO AMAZONAS -0.7568694 -0.4650969 -1.8614255 -0.3405284 -1.4160776 -0.5268008
## AMBO           1.9002588  0.4136115  0.1684413 -0.4813803  0.1304697 -0.3309393
##                    Axis7       Axis8       Axis9      Axis10      Axis11
## ABANCAY       -0.5672185 -0.49320512 -0.32000416 -0.62199816 -0.17433866
## ACOBAMBA       0.6262939  0.20368922 -0.07566004  0.10789738  0.55611620
## ACOMAYO        0.3558107 -0.02599267 -0.69317540 -0.12428901  0.21103518
## AIJA          -1.1117521 -0.04258963 -0.21010040 -0.03939282  0.12359601
## ALTO AMAZONAS  0.1962207 -0.55845951  0.22637417  0.11677574 -0.08454269
## AMBO           0.0106995 -0.19401066  0.32137963 -0.06387590  0.29155440
##                     Axis12      Axis13
## ABANCAY        0.300250564  0.16745855
## ACOBAMBA      -0.051762227 -0.02850826
## ACOMAYO        0.003353471 -0.07537753
## AIJA          -0.097778789 -0.19194143
## ALTO AMAZONAS  0.154516096 -0.22993921
## AMBO          -0.017675385  0.04726788
head(PCA_1ola$li)
##                    Axis1      Axis2      Axis3      Axis4      Axis5      Axis6
## ABANCAY       -3.4919948  2.7915258 -3.1313116  0.2103763  0.7547921 -0.4381602
## ACOBAMBA       2.0963553 -0.4770106 -0.1553921  2.2149076  0.3234739  1.6660686
## ACOMAYO        2.1663012  0.7270715  0.3506939 -0.1628312  0.2269999  0.2316588
## AIJA           1.8371843  0.4596678  0.6078255 -0.1877961  0.4789561  0.3981455
## ALTO AMAZONAS -0.7568694 -0.4650969 -1.8614255 -0.3405284 -1.4160776 -0.5268008
## AMBO           1.9002588  0.4136115  0.1684413 -0.4813803  0.1304697 -0.3309393
##                    Axis7       Axis8       Axis9      Axis10      Axis11
## ABANCAY       -0.5672185 -0.49320512 -0.32000416 -0.62199816 -0.17433866
## ACOBAMBA       0.6262939  0.20368922 -0.07566004  0.10789738  0.55611620
## ACOMAYO        0.3558107 -0.02599267 -0.69317540 -0.12428901  0.21103518
## AIJA          -1.1117521 -0.04258963 -0.21010040 -0.03939282  0.12359601
## ALTO AMAZONAS  0.1962207 -0.55845951  0.22637417  0.11677574 -0.08454269
## AMBO           0.0106995 -0.19401066  0.32137963 -0.06387590  0.29155440
##                     Axis12      Axis13
## ABANCAY        0.300250564  0.16745855
## ACOBAMBA      -0.051762227 -0.02850826
## ACOMAYO        0.003353471 -0.07537753
## AIJA          -0.097778789 -0.19194143
## ALTO AMAZONAS  0.154516096 -0.22993921
## AMBO          -0.017675385  0.04726788
dim(PCA_1ola$li)
## [1] 196  13

ENTONCES COMO DEFINIMOS QUE NOS QUEDAMOS LAS 2 PRIMERAS COMPONENTES, SELECCIONAMOS SOLO DICHAS COMPONENTES.

output_scores_1ola <-
  as.tibble(PCA_1ola$li) %>% 
  dplyr::select(sprintf("Axis%1$s",1:2))

View(output_scores_1ola)

CALCULAR UN INDICE A PARTIR DE COMPONENTES PRINCIPALES:

I = ((PCA_1ola$li$Axis1*PCA_1ola$eig[1])+(PCA_1ola$li$Axis2*PCA_1ola$eig[2]))/(PCA_1ola$eig[1]+PCA_1ola$eig[2])

Indice_scores_1ola <-dplyr:: bind_cols(output_scores_1ola, Indice = I) %>% 
  cbind(Data_1_Ola["PROVINCIA"])

View(Indice_scores_1ola)
  • pca\(li\)Axis1: Representa los scores o puntuaciones de la primer componente
  • pca$eig[1]: son los autovalores de la primera componente

EL CALCULO DEL INDICE SIRVE PARA EL MONITOREO, POR EJEMPLO PODEMOS MONITOREAR LA OCURRENCIA DE UN FENÓMENO A PARTIR DE UN SOLO INDICADOR #ESTE INDICE SE PUEDE ESTANDARIZAR CON 0 Y 1, DONDE 0 TE INDICA VALORES BAJO DE PROBABILIDAD DE INCENDIOS POR EJEMPLO Y 1 REPRESENTA VALORES ALTOS #DE INCENDIOS

ANÁLISIS CLUSTER:

ANÁLISIS CLUSTER AGLOMERATIVO:

CÁLCULO DE LAS DISTANCIAS: Para realizar un análisis cluster aglomerativo, se debe primero determinar las distancias en función al método, en este caso estamos utilizando las distancias euclidianas

distancia_1ola <- dist(output_scores_1ola, method = "euclidean") #Colocamos los scores del PCA resultante

Para calcular la dimension de la matriz de distancia es de la siguiente manera:

as.matrix(distancia_1ola) %>% dim()
## [1] 196 196

OBTENIENDO EL DENDOGRAMA:

Para calcular el dendograma se aplicará la función “hclust”, en donde se aplicará le método de agrupamiento “ward.D2”, dado que las observaciones de las variables se encuentran muy juntos:

hmodel_1ola <- hclust(distancia_1ola, method ="ward.D2")

PLOTEANDO EL DENDOGRAMA: El ploteo del dendograma se realizará con R base:

par("bg") 
## [1] "white"
par(bg = "black") #Color del fondo
plot(hmodel_1ola,
     main="Dendograma",# Título
     las = 1,                 # Rotar las etiquetas de los ejes de manera horizontal
     xlab = "Cluster",          # Etiqueta del eje x
     ylab = "Alturas",          # Etiqueta del eje y
     cex.main = 1.5,          # Tamaño del título
     cex.sub = 1.2,           # Tamaño del subtítulo
     cex.lab = 1.2,           # Tamaño de las etiquetas de los ejes X e Y
     cex.axis = 1,            # Tamaño de las etiquetas de los ticks de los ejes
     font.sub = 1,            # Estilo de fuente del título sin formato
     font.main  = 2,          # Estilo de fuente del subtítulo en negrita
     font.axis = 3,           # Estilo de fuente de los ejes X e Y en cursiva
     font.lab  = 4,           # Estilo de fuente de los ticks de los ejes en negrita y cursiva
     pch  = 21,               # Estilo del símbolo de los puntos del grafico
     bg = "red",              # Color de fondo del simbolo
     col = "white",           # Color del borde del simbolo
     cex = 1.2,               # Tamaño del símbolo
     lwd = 1.2,               # Ancho del borde
     col.main = "Yellow",      # Color del título
     col.sub = "blue",        # Color del subtítulo
     col.lab = "orange",      # Color de las etiqetas de los ejes
     col.axis = "maroon4",    # Color de las etiquetas de los ticks
     fg = "green") 

El eje de las alturas es un valor creciente, es decir aumenta en función a cómo se va construyendo los cluster, debido a que la cohesion disminuye cada vez que aglomeramos los cluster o exigimos que se unan.

PROCESO DE AGRUPAMIENTO INDICANDO DISTANCIAS:

Para obtener las alturas aplicaremos el siguiente código:

head(hmodel_1ola$height)
## [1] 0.02207779 0.04089096 0.04127458 0.04171793 0.04483989 0.04625527

Ploteando el gráfico de las alturas:

plot.new()

rect(par("usr")[1], par("usr")[3],
     par("usr")[2], par("usr")[4],
     col = "black") 

par(new = TRUE)

plot(hmodel_1ola$height,
     main="Gráfico de alturas",# Título
     las = 1,                 # Rotar las etiquetas de los ejes de manera horizontal
     xlab = "Index",          # Etiqueta del eje x
     ylab = "Alturas",          # Etiqueta del eje y
     cex.main = 1.5,          # Tamaño del título
     cex.sub = 1.2,           # Tamaño del subtítulo
     cex.lab = 1.2,           # Tamaño de las etiquetas de los ejes X e Y
     cex.axis = 1,            # Tamaño de las etiquetas de los ticks de los ejes
     font.sub = 1,            # Estilo de fuente del título sin formato
     font.main  = 2,          # Estilo de fuente del subtítulo en negrita
     font.axis = 3,           # Estilo de fuente de los ejes X e Y en cursiva
     font.lab  = 4,           # Estilo de fuente de los ticks de los ejes en negrita y cursiva
     pch  = 21,               # Estilo del símbolo de los puntos del grafico
     bg = "red",              # Color de fondo del simbolo
     col = "black",           # Color del borde del simbolo
     cex = 1.2,               # Tamaño del símbolo
     lwd = 1.2,               # Ancho del borde
     col.main = "Orange",      # Color del título
     col.sub = "blue",        # Color del subtítulo
     col.lab = "orange",      # Color de las etiqetas de los ejes
     col.axis = "maroon4",    # Color de las etiquetas de los ticks
     fg = "green")            # Color de la caja
lines(hmodel_1ola$height, col = "yellow") 

#Para agregar una grilla es con la función "Grid"

grid(nx = NULL, ny = NULL,    # Agrega una grilla a los ejes X e Y
     lty = 2,                 # Tipo de línea
     col = "gray",            # Color
     lwd = 1)                 # Ancho de línea

par(new = TRUE)

Este proceso de inclusión de clúster según el ploteo resultante es hasta llegar a un punto en donde los clúster se han unido coherentemente, pero luego le ha costado fundirse con otro clúster que es bastante diferente, entonces observamos un quiebre. El objetivo es encontrar ese punto de quiebre para cortar y a partir de ahí, contamos cuántos cluster tenemos,en este caso vemos que se puede observar un punto de quiebre pronunciado en la observación 194, pero si nos fijamos mása detalle, se nota un quiebre de inicio en la observación 191.

Bajo esa índole analizaremos qué valor de altura le corresponde a la observación 191:

hmodel_1ola$height[191]
## [1] 10.62108

Entonces para una altura de 10.621 se obtienen 6 clúster.

Imagen

Entonces ploteamos el dendograma resultante

par("bg")
## [1] "white"
par(bg = "black")
plot(hmodel_1ola,
     main="Dendograma",# Título
     las = 1,                 # Rotar las etiquetas de los ejes de manera horizontal
     xlab = "Cluster",          # Etiqueta del eje x
     ylab = "Alturas",          # Etiqueta del eje y
     cex.main = 1.5,          # Tamaño del título
     cex.sub = 1.2,           # Tamaño del subtítulo
     cex.lab = 1.2,           # Tamaño de las etiquetas de los ejes X e Y
     cex.axis = 1,            # Tamaño de las etiquetas de los ticks de los ejes
     font.sub = 1,            # Estilo de fuente del título sin formato
     font.main  = 2,          # Estilo de fuente del subtítulo en negrita
     font.axis = 3,           # Estilo de fuente de los ejes X e Y en cursiva
     font.lab  = 4,           # Estilo de fuente de los ticks de los ejes en negrita y cursiva
     pch  = 21,               # Estilo del símbolo de los puntos del grafico
     bg = "red",              # Color de fondo del simbolo
     col = "white",           # Color del borde del simbolo
     cex = 1.2,               # Tamaño del símbolo
     lwd = 1.2,               # Ancho del borde
     col.main = "Yellow",      # Color del título
     col.sub = "blue",        # Color del subtítulo
     col.lab = "orange",      # Color de las etiqetas de los ejes
     col.axis = "maroon4",    # Color de las etiquetas de los ticks
     fg = "green") 
rect.hclust(hmodel_1ola, k = 6,
            border = 3:4)

Entonces como resultado tenemos que el K óptimo es 6.

ANÁLISIS CLUSTER DE PARTICIONAMIENTO:

Para el análisis cluster de particionamiento se aplicará el método de Kmeans++, pero antes de eso evaluaremos otra forma de determinar el K óptimo.

Existen 3 formas principales para realizar el cálculo del K óptimo, estos son los siguientes:

a) MÉTODO DE LA SILUETA:

Este gráfico te indica cuánto es el valor óptimo, en función al que presenta mayor valor de silueta:

factoextra::fviz_nbclust( 
  output_scores_1ola, kmeans, method = "silhouette",#Usaremos el metodo de  la silueta
  k.max = 20, #Establecemos el número de cluster máximo.
  linecolor = "blue",
  print.summary = TRUE
)

Según este gráfico, el que presenta mayor valor de silueta es el número de clúster 2, es decir el K óptimos según este método sería de 2.

b) MÉTODO DEL CODO:

Este gráfico te indica cuánto es el valor óptimo, en función al punto de inflexión:

factoextra::fviz_nbclust( 
  output_scores_1ola, kmeans, method = "wss",#Usaremos el metodo de  sumas de cuadrados
  k.max = 20 #Establecemos el número de cluster máximo.
)

Según este gráfico, vemos que el punto de inflexión se encuentra en el cluster 2, por lo que el K óptimo sería 2.

c) MÉTODO DE LOS INDICADORES:

Este método te determina un conjunto de indicadores que te comentan la cantidad óptima de clusteres a formar:

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 8 proposed 2 as the best number of clusters 
## * 5 proposed 3 as the best number of clusters 
## * 2 proposed 4 as the best number of clusters 
## * 1 proposed 8 as the best number of clusters 
## * 3 proposed 9 as the best number of clusters 
## * 1 proposed 15 as the best number of clusters 
## * 3 proposed 20 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  2 
##  
##  
## *******************************************************************
## Among all indices: 
## ===================
## * 2 proposed  0 as the best number of clusters
## * 1 proposed  1 as the best number of clusters
## * 8 proposed  2 as the best number of clusters
## * 5 proposed  3 as the best number of clusters
## * 2 proposed  4 as the best number of clusters
## * 1 proposed  8 as the best number of clusters
## * 3 proposed  9 as the best number of clusters
## * 1 proposed  15 as the best number of clusters
## * 3 proposed  20 as the best number of clusters
## 
## Conclusion
## =========================
## * According to the majority rule, the best number of clusters is  2 .

Este gráfico te indica que existe 8 indicadores que determinan que el K óptimo sea 2.

ANALISIS DE CLUSTER CON PCA

Vamos a generar el gráfico entre las dos primeras componentes, para esta ocasión se colocará la data estandarizada y no los scores resultantes, dado que el modelo corre por debajo el PCA:

#factoextra::fviz_cluster(
  #object = model_1ola , #Insertamos la salida del modelo
  #data = scale_data_1ola, #Insertamos la data original estandarizada
  #ellipse.type = "euclid", #colocamos el tipo de elipse (pUEDE SER TAMBIEN CONVEX)
  #outlier.color = "black" #colorea aquellas observaciones alejadas
#)

Hay que resaltar que este gráfico no es un método de clusterización, lo que estamos observando es cuán heterogéneo son los individuos para las 2 primeras componentes.

APLICACIÓN DEL MÉTODO DE PARTICIONAMIENTO KMEANS++:

Para realizar el método de particionamiento Kmeans ++, tenemos que establecer una semilla, dado que el centroide central es aleatorio. En este caso, comenzaremos con un K óptimo de 2:

set.seed(2021) 
model1_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola, #Colocar nuestros datos estandarizados o nuestros scores de las componentes que hemos retenido
    k = 2, #Establecemos el k óptimo
    start = "random", #Start es diferente que nstart, esto significa que te da a escoger que el primer centroide sea seleccionado de manera ramdom o que el centroide tenga un lugar especifico.
    iter.max = 100,#Se coloca el número de iteraciones
    nstart = 200,# Todos los modelos que se corran paralelamente, van a tener los mismo resultados. 
    algorithm = "Hartigan-Wong", #Escogemos el algoritmo de kmeans.
    trace = F  #Es como una especie de reporte que te genera acerca de las iteraciones, en este caso ponemos FALSE. 
  )

Adicionalmente podemos tener estas informaciones:

model1_1ola$withinss #Te muestra la suma de cuadrados intracuster
## [1] 246.2101 448.0365
model1_1ola$tot.withinss #Te muestra la suma de cuadrados intracluster total
## [1] 694.2466
model1_1ola$betweenss #Te muestra la suma de cuadrados entre cluster (intercluster)
## [1] 890.6565
model1_1ola$totss #Te muestra la suma de cuadrados entre cluster total
## [1] 1584.903
model1_1ola$cluster #Te muestra los cluster resultantes
##   [1] 2 1 1 1 2 1 1 1 1 1 1 2 2 1 1 1 1 1 2 2 1 1 1 1 1 2 1 1 2 2 1 1 1 1 2 2 1
##  [38] 2 1 1 2 1 1 2 1 2 2 2 2 2 1 1 1 1 1 1 1 1 1 2 1 2 1 1 2 1 1 1 1 1 1 1 2 1
##  [75] 1 1 1 1 2 1 2 1 1 1 2 2 1 2 2 2 1 2 1 2 1 1 2 2 2 2 2 2 1 1 1 1 1 1 2 1 1
## [112] 2 2 1 1 1 2 1 2 1 2 1 2 1 1 2 2 2 1 1 1 1 2 1 2 2 1 2 1 1 2 1 1 1 1 2 2 1
## [149] 1 2 2 2 1 1 1 1 1 1 1 1 2 1 1 2 1 1 2 1 1 1 2 1 1 2 2 2 2 2 1 1 1 2 2 2 1
## [186] 1 2 1 1 1 1 2 1 1 1 2
model1_1ola$centers #Te muestra los centroides en coordenadas multidimensionales, estos centroides son los valores promedios de cada variable u observacion, son los centrioides finales.
##       Axis1       Axis2
## 1  1.641601  0.04570870
## 2 -2.765986 -0.07701602
model1_1ola$inicial.centers #Te muestra los centroides iniciales
##         Axis1     Axis2
## 135 -1.244399 -1.537883
## 65  -3.062261  1.452068
model1_1ola$size #Te muestra la cantidad de individuos por cada cluster, en total debe salir 196
## [1] 123  73

Por consiguiente, obtendremos el gráfico de silueta:

Valor de la silueta

sil_1ola <-cluster::silhouette(##Te calcula el valor de la silueta para cada individuo
  model1_1ola$cluster, #Insertar los cluster del modelo
  dist(output_scores_1ola)) #Insertar la matriz de distancia.

Gráfico de la silueta

F1<-factoextra::fviz_silhouette(sil_1ola) +
  coord_flip ()+
  theme_bw()
##   cluster size ave.sil.width
## 1       1  123          0.63
## 2       2   73          0.35
plot(F1)

En el eje “x” aparecen las variables, cada cluster tiene su valor de silueta promedio. Podemos decir que el valor promedio de silueta es de 0.52, lo que significa que presenta un alto valor, por lo que el método se ha clusterizado de manera correcta. Así mismo vemos que presenta valores de silueta negativos para algunas observaciones.

Entonces otro método que te determinará el K óptimo es mediante los gráficos de la silueta, donde aquel resultado que presente la menor cantidad de observaciones negativas, será el resultado del K óptimo, para eso evaluaremos los siguientes clústeres:

Con 3 clúster

set.seed(2021) 
model2_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola, 
    k = 3, 
    start = "random", 
    iter.max = 100,
    nstart = 200, 
    algorithm = "Hartigan-Wong", 
    trace = F  
  )


#VALORES DE SILUETA:
sil2_1ola <-cluster::silhouette(
  model2_1ola$cluster, 
  dist(output_scores_1ola)) 


#GRAFICANDO LA SILUETA:

F2<-factoextra::fviz_silhouette(sil2_1ola) +
  coord_flip () +
  theme_bw()
##   cluster size ave.sil.width
## 1       1   54          0.28
## 2       2   33          0.37
## 3       3  109          0.61

Con 4 clúster

set.seed(2021)
model3_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola, 
    k = 4, 
    start = "random", 
    iter.max = 100,
    nstart = 200,
    algorithm = "Hartigan-Wong", 
    trace = F  
  )


#VALORES DE SILUETA:
sil3_1ola <-cluster::silhouette(
  model3_1ola$cluster, 
  dist(output_scores_1ola)) 


#GRAFICANDO LA SILUETA:

F3<-factoextra::fviz_silhouette(sil3_1ola) +
  coord_flip () +
  theme_bw()
##   cluster size ave.sil.width
## 1       1  104          0.62
## 2       2   21          0.23
## 3       3   25          0.25
## 4       4   46          0.37

Con 5 clúster

set.seed(2021) 

model4_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola, 
    k = 5, 
    start = "random", 
    iter.max = 100,
    nstart = 200,
    algorithm = "Hartigan-Wong", 
    trace = F  
  )


#VALORES DE SILUETA:
sil4_1ola <-cluster::silhouette(
  model4_1ola$cluster, 
  dist(output_scores_1ola))  


#GRAFICANDO LA SILUETA:

F4<-factoextra::fviz_silhouette(sil4_1ola) +
  coord_flip () +
  theme_bw()
##   cluster size ave.sil.width
## 1       1   98          0.58
## 2       2    3          0.23
## 3       3   48          0.33
## 4       4   35          0.36
## 5       5   12          0.29

Con 6 clúster

set.seed(2021) 
model5_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola,  
    k = 6, 
    start = "random", 
    iter.max = 100,
    nstart = 200,
    algorithm = "Hartigan-Wong", 
    trace = F  
  )


#VALORES DE SILUETA:
sil5_1ola <-cluster::silhouette(
  model5_1ola$cluster,
  dist(output_scores_1ola)) 



#GRAFICANDO LA SILUETA:

F5<-factoextra::fviz_silhouette(sil5_1ola) +
  coord_flip () +
  theme_bw()
##   cluster size ave.sil.width
## 1       1   17          0.30
## 2       2   25          0.45
## 3       3   95          0.56
## 4       4    3          0.21
## 5       5   19          0.30
## 6       6   37          0.37

Con 7 clúster

set.seed(2021) 

model6_1ola <- 
  LICORS::kmeanspp(
    data = output_scores_1ola, 
    k = 7, 
    start = "random", 
    iter.max = 100,
    nstart = 200,
    algorithm = "Hartigan-Wong", 
    trace = F  
  )


#VALORES DE SILUETA:
sil6_1ola <-cluster::silhouette(
  model6_1ola$cluster, 
  dist(output_scores_1ola)) 


#GRAFICANDO LA SILUETA:

F6<-factoextra::fviz_silhouette(sil6_1ola) +
  coord_flip () +
  theme_bw()
##   cluster size ave.sil.width
## 1       1   18          0.25
## 2       2   46          0.33
## 3       3    3          0.21
## 4       4   59          0.39
## 5       5   29          0.39
## 6       6   16          0.29
## 7       7   25          0.43

Por consiguiente, plotearemos en un solo gráficos las siluetas resultantes

require(ggpubr)
## Loading required package: ggpubr
ggpubr::ggarrange(F1,F2,F3,F4,F5,F6)

Observamos que el gráfico que presenta menores valores de siluetas negativas el aquel que presenta 6 clústeres, por lo que el K óptimo es de 6.

ANÁLISIS DE CLUSTER CON PCA:

factoextra::fviz_cluster(
  object = model5_1ola , #Insertamos la salida del modelo
  data = scale_data_1ola, #Insertamos la data original estandarizada
  ellipse.type = "convex", #colocamos el tipo de figura
  outlier.color = "red" ,#colorea aquellas observaciones alejadas
  ggtheme = theme_gray(),
  main = "                    
  Ploteo de los clústeres")

En este gráfico podemos observar el análisis de los 2 componentes con los 6 clústeres formados, hay que resaltar que el objetivo de este gráfico es observar cuán heterogéneo son los individuos para las 2 primeras componentes, así mismo hay que recordar que se coloca la data estandarizada y no los scores, debido a que este modelo te genera por debajo el PCA.

Resultado del análisis clúster

En la parte final, vamos a obtener una tabla resultante:

Output2_1ola<- cbind(Data_1_Ola["PROVINCIA"]) %>% 
  dplyr::bind_cols(output_scores_1ola, cluster = model5_1ola$cluster) %>% 
  cbind(Indice_scores_1ola["Indice"]) %>% 
  cbind(Data_1_Ola) %>% 
  dplyr::select(c( ,-6))

head(Output2_1ola)
##                   PROVINCIA      Axis1      Axis2 cluster     Indice
## ABANCAY             ABANCAY -3.4919948  2.7915258       2 -2.1604852
## ACOBAMBA           ACOBAMBA  2.0963553 -0.4770106       3  1.5510461
## ACOMAYO             ACOMAYO  2.1663012  0.7270715       3  1.8613212
## AIJA                   AIJA  1.8371843  0.4596678       3  1.5452816
## ALTO AMAZONAS ALTO AMAZONAS -0.7568694 -0.4650969       6 -0.6950414
## AMBO                   AMBO  1.9002588  0.4136115       3  1.5852308
##               ElevacionPromedio IDH_2019 FallMasc_1ola_mill_hab
## ABANCAY                    3100     0.52               749.3131
## ACOBAMBA                   3300     0.31               370.9153
## ACOMAYO                    3800     0.30               353.3694
## AIJA                       2800     0.37               252.0479
## ALTO AMAZONAS              1200     0.44              1211.2267
## AMBO                       3300     0.38               323.0413
##               FallFem_1ola_mill_hab CasosMasc_1ola_mill_hab
## ABANCAY                    527.2944               11702.236
## ACOBAMBA                   123.6384                3090.961
## ACOMAYO                    353.3694                2120.216
## AIJA                         0.0000                5923.125
## ALTO AMAZONAS              754.0087                8037.412
## AMBO                       255.0326                3111.398
##               CasosFem_1ola_mill_hab Adult_15_64_mill_hab Mayores_65_mill_hab
## ABANCAY                    13774.410             691505.0            85662.22
## ACOBAMBA                    4030.613             303371.6            44262.56
## ACOMAYO                     1908.195             512915.7            96929.22
## AIJA                        3780.718             491115.3           116824.20
## ALTO AMAZONAS               6072.176             616666.8            54577.39
## AMBO                        3621.464             552961.8            96079.30
##               Pobreza_mill_hab PobrezaExt_mill_hab UCI_1ola_cant_hosp
## ABANCAY               265579.3            84670.16           8.666667
## ACOBAMBA              247370.7            94642.38           0.000000
## ACOMAYO               386451.8           134761.65           0.000000
## AIJA                  323041.0           130509.14           0.000000
## ALTO AMAZONAS         441580.2           253656.06           3.000000
## AMBO                  357757.2           239162.64           0.000000
##               UCIN_1ola_cant_hosp VENT_1ola_cant_hosp
## ABANCAY                      78.0            8.666667
## ACOBAMBA                      6.0            0.000000
## ACOMAYO                       0.0            0.000000
## AIJA                          0.0            0.000000
## ALTO AMAZONAS                33.5            3.000000
## AMBO                          0.0            0.000000

Ploteando el resultado

Una vez obtenida la tabla resultante, se procederá a crear un mapa con el resultado de la clasterización:

#Cargamos el Shp de Provincias:
SHP_1ola <- st_read("Provincia/PROVINCIA_FINAL.shp")
## Reading layer `PROVINCIA_FINAL' from data source `C:\Users\USER\Desktop\clase\NOVENO CICLO\PROGRAMACION II\ARCHIVOS ULTIMO TRABAJO DE R\Provincia\PROVINCIA_FINAL.shp' using driver `ESRI Shapefile'
## Simple feature collection with 196 features and 1 field
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -81.32823 ymin: -18.35093 xmax: -68.65228 ymax: -0.03860597
## geographic CRS: WGS 84
#Unimos la tabla resultante con el shp de provincia:
SHP_FINAL_1OLA<-Output2_1ola %>% 
  left_join(SHP_1ola, by = "PROVINCIA")


#Ploteando

Mapa_1ola<-st_sf(SHP_FINAL_1OLA) %>% 
  dplyr::select(cluster,PROVINCIA) 

ggplot(Mapa_1ola) +
  geom_sf(aes(fill = cluster))+
  labs(title = "Mapa de clusterización de la primera ola",
       caption = "Elaboración propia",
       x="Longitud",
       y="Latitud") +
  scale_fill_gradientn("Clúster", colours = c("orange","yellow","green","blue","Turquoise","pink"))+
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 9, face = "bold"))

VALIDACIÓN DEL ANÁLISIS CLUSTER

INDICE DE DAVIES-BOULDING

Este índice te dice que mientras sea más cercano a cero, el proceso de clusterización sería óptimo.

grupo_1ola <- model5_1ola$cluster

index <- clusterSim::index.DB(
  output_scores_1ola, #Ingresamos las datas estandarizadas o los scores retenidos
  grupo_1ola, #Ingresamos los cluster resultantes
  centrotypes = "centroids" # ingresamos el centro que se medira, en este caso de los centroides
)


#OBTENIENDO EL INDICE:
index$DB
## [1] 0.8835893

En este caso vemos que el valor supera al 1, pero hay que recordar que solo es un indicador que no debe sesgar el resultado, por lo que validamos el proceso de clusterización mediante el gráfico de la silueta.

INDICE DE DUNN

Este índice te dice que mientras el valor sea más grande, el proceso de clusterización sería óptimo.

clValid::dunn(
  Data = output_scores_1ola, #Ingresamos las datas estandarizadas o los scores retenidos
  cluster = grupo_1ola, #Ingresamos los cluster resultantes
  distance = NULL 
)
## [1] 0.0428975

El índice de DAM tiene que ser lo más grande posible, es al contrario del anterior índice de DAVIES-BOULDING, en este caso vemos que el valor es cercano a 0, pero de la misma manera, solo es un indicador que no debe sesgar el resultado, por lo que validamos el proceso de clusterización mediante el gráfico de la silueta.

Análisis cluster

Gráficos

Cargar las variables almacenadas en el archivo Rdata

load("Analisis2.RData")

Fallecidos masculinos:

G1

Altitudes:

G2

Fallecidos masculinos:

G3

Fallecidos masculinos:

G4

Fallecidos masculinos:

G5

Fallecidos masculinos:

G6

Población con edad entre 15 a 64 años:

G7

Pobreza:

G8

Pobreza extrema:

G9

Población con edad superior a 65 años:

G10

Unidad de cuidados intensivos(UCI), por hospital:

G11

Unidad de cuidados intensivos(UCI) intermedio, por hospital:

G12

Cantidad de ventiladores por hospital:

G13